#include "base.h"

int main(int argc,char **argv)
{


          //////////////////
          // Input   ///////
          /////////////////

	  int NOa2;
	  int Nf;
          double beta, Uplaq;
          char dir[1024];
          char gauge_fixing[1024];
          char gauge_action[1024];
          char momentum_order[1024];
	  int beta100;
	  int Nr=2;
	  int Nmass;

          ////// Parsing command line

          for(int i=1; i<argc; i++) {
            if(strncmp(argv[i],"-Nclust",7)==0)  sscanf(argv[i],"-Nclust=%d",&Nclust);
            if(strncmp(argv[i],"-Nf",3)==0)  sscanf(argv[i],"-Nf=%d",&Nf);
            if(strncmp(argv[i],"-dir",4)==0)  sscanf(argv[i],"-dir=%s",dir);
            if(strncmp(argv[i],"-gauge_fixing",13)==0)  sscanf(argv[i],"-gauge_fixing=%s",gauge_fixing);
            if(strncmp(argv[i],"-gauge_action",13)==0)  sscanf(argv[i],"-gauge_action=%s",gauge_action);
            if(strncmp(argv[i],"-T",2)==0)  sscanf(argv[i],"-T=%d",&T);
            if(strncmp(argv[i],"-L",2)==0)  sscanf(argv[i],"-L=%d",&L);
            if(strncmp(argv[i],"-beta",5)==0)  sscanf(argv[i],"-beta=%lf",&beta);
            if(strncmp(argv[i],"-Uplaq",6)==0)  sscanf(argv[i],"-Uplaq=%lf",&Uplaq);
            if(strncmp(argv[i],"-NOa2",5)==0)  sscanf(argv[i],"-NOa2=%d",&NOa2);
            if(strncmp(argv[i],"-momentum_order",15)==0)  sscanf(argv[i],"-momentum_order=%s",momentum_order);
            if(strncmp(argv[i],"-Nmass",6)==0)  sscanf(argv[i],"-Nmass=%d",&Nmass);
            if(strncmp(argv[i],"-beta100",8)==0)  sscanf(argv[i],"-beta100=%d",&beta100);

          }

   /////////////////////////
  //// Momenta set
  //////////////////////////////

  char momentafile[1024]; sprintf(momentafile,"../INPUT/mom_%d.dat",beta100);
  read_momenta_data(momentafile);
  double *nn=(double*)malloc(sizeof(double)*(Nmom));
  compute_nn(momentum_order,nn);
  Momentum(momentum_order);
  // equivalent momenta
  Nmom_eq=Number_of_Equivalent_Momenta(SinP2);
  Equivalent_Momenta(SinP2);

	printf(" Nmom %d NmomEq %d\n",Nmom, Nmom_eq);
  // IO

   char file[1024];
   char Oa2print[1024];
   char dir_qav[1024];

    //////////////////////////////////
    // Read Zq 
     /////////////////////////////
     printf("====================\n");
     printf(" READ   \n");

    // Note that the data is complex
    // But from now on I'm only interested in the real part

        double ****Zq=malloc(sizeof(double***)*(Nclust+1));
        for (int iclust=0; iclust<Nclust+1; iclust++) Zq[iclust]=malloc(sizeof(double**)*(Nmass));
        for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) Zq[iclust][imass]=malloc(sizeof(double*)*(Nmom));
        for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) Zq[iclust][imass][imom]=malloc(sizeof(double)*(Nr));


    complex *aux_clust=malloc(sizeof(complex)*(Nclust+1));

    for (int ir=0; ir<2; ir++) for (int imass=0; imass<Nmass; imass++)
    {
    sprintf(file,"Zq_r%d_m%d",ir,imass);
    for (int imom=0; imom<Nmom; imom++) {
	    read_input_onemomenta_complex(dir,file,aux_clust,imom);
	    for (int iclust=0; iclust<Nclust+1; iclust++) Zq[iclust][imass][imom][ir]=creal(aux_clust[iclust]);
     }
     }

     printf("====================\n");
     printf(" PROMEDIUM r0 r1   \n");
     double ****Zqr=malloc(sizeof(double***)*(Nclust+1));
     for (int iclust=0; iclust<Nclust+1; iclust++) Zqr[iclust]=malloc(sizeof(double**)*(Nmass));
     for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) Zqr[iclust][imass]=malloc(sizeof(double*)*(Nmom));
      for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) Zqr[iclust][imass][imom]=malloc(sizeof(double)*(NOa2));

   for (int iclust=0; iclust<Nclust+1; iclust++) for (int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++)
        {
                Zqr[iclust][imass][imom][0]=(Zq[iclust][imass][imom][0]+Zq[iclust][imass][imom][1])/2.;
        }


    printf("====================================================\n");
     printf("O(a2g2) corrections \n ");


     //////////////// substracting terms     
     double ***CorOa2_Zq=malloc(sizeof(double**)*(NOa2));
     for (int iOa2=0; iOa2<NOa2; iOa2++) CorOa2_Zq[iOa2]=malloc(sizeof(double*)*(Nmom));
     for (int iOa2=0; iOa2<NOa2; iOa2++)  for (int imom=0; imom<Nmom; imom++) CorOa2_Zq[iOa2][imom]=malloc(sizeof(double)*(1));
     double ***CorOa2_Lambda=malloc(sizeof(double**)*(NOa2));
     for (int iOa2=0; iOa2<NOa2; iOa2++) CorOa2_Lambda[iOa2]=malloc(sizeof(double*)*(Nmom));
     for (int iOa2=0; iOa2<NOa2; iOa2++)  for (int imom=0; imom<Nmom; imom++) CorOa2_Lambda[iOa2][imom]=malloc(sizeof(double)*(5));

     //iOa2=0 -> uncorrected
     for (int imom=0; imom<Nmom; imom++){
     CorOa2_Zq[0][imom][0]=0.;
     for (int igamma=0; igamma<5 ; igamma++) CorOa2_Lambda[0][imom][igamma]=0.;
     }
     //iOa2=1 corrected
     //iOa2=2 corrected boosted
     double g2;
     for (int iOa2=1; iOa2<NOa2; iOa2++)
     {
        for (int imom=0; imom<Nmom; imom++)
        {
          if(iOa2==1) g2=6.0/beta;
          if(iOa2==2) g2=(6.0/beta)/Uplaq;
          correction_Oa2g2(CorOa2_Zq[iOa2][imom],CorOa2_Lambda[iOa2][imom],g2,SinP2[imom],SinP4[imom],nn[imom],gauge_action,gauge_fixing);
        }

     }
      ////////////////// Zq corrected
      for (int iclust=0; iclust<Nclust+1; iclust++) for (int imom=0; imom<Nmom; imom++)
      for (int imass=0; imass<Nmass; imass++)
      {
         for (int iOa2=1; iOa2<NOa2; iOa2++) Zqr[iclust][imass][imom][iOa2]=Zqr[iclust][imass][imom][0]-CorOa2_Zq[iOa2][imom][0];
      }

     //////////////////////////////////////////////
     //     EQUIVALENT MOMENTA
     ///////////////////////////////////////////////
     printf("==================================\n");
     printf(" PROMEDIUM OVER EQUIVALENT SINP2 \n");

     double ****Zqr_eq=malloc(sizeof(double***)*(Nclust+1));
     for (int iclust=0; iclust<Nclust+1; iclust++) Zqr_eq[iclust]=malloc(sizeof(double**)*(Nmass));
     for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) Zqr_eq[iclust][imass]=malloc(sizeof(double*)*(Nmom_eq));
     for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom_eq; imom++)  Zqr_eq[iclust][imass][imom]=malloc(sizeof(complex)*(NOa2));

     double *aux_mom=malloc(sizeof(double)*Nmom);
     double *aux_eq_mom=malloc(sizeof(double)*Nmom_eq);

   for (int iOa2=0; iOa2<NOa2; iOa2++)
   for (int imass=0; imass<Nmass; imass++)
    for (int iclust=0; iclust<Nclust+1; iclust++)
   {
    for (int imom=0; imom<Nmom; imom++) aux_mom[imom]=Zqr[iclust][imass][imom][iOa2];
    Average_over_Equivalent_Momenta(SinP2,aux_mom,aux_eq_mom);
    for (int imom=0; imom<Nmom_eq; imom++)  Zqr_eq[iclust][imass][imom][iOa2]=aux_eq_mom[imom];

   }

   ///////////
   // write output
   ///////////////
  printf("=========================\n");
  printf(" WRITING OUTPUT\n");
  
   sprintf(dir_qav,"%s/QAV/",dir);
   double *aux_clust_real=malloc(sizeof(double)*(Nclust+1));

        // Zq RPROM QAVERAGE
    for (int iOa2=0; iOa2<NOa2; iOa2++)
    for (int imass=0; imass<Nmass; imass++)
     {
        whichOa2(iOa2, Oa2print);
        sprintf(file,"Zq_%s_m%d",Oa2print,imass);
        printf("Writting %s file output in %s \n",file,dir_qav);

         for (int imom=0; imom<Nmom_eq; imom++){
          for (int iclust=0; iclust<Nclust+1; iclust++)
                aux_clust_real[iclust]=Zqr_eq[iclust][imass][imom][iOa2];
          write_output_onemomenta_real(dir_qav,file,aux_clust_real,SinP2_eq,imom);
         }

         sprintf(file,"%s/Z_ALLP/Zq_%s_m%d.dat",dir_qav,Oa2print,imass);
         FILE *fout=fopen(file,"w");
         if(fout == (FILE *)NULL) {
              fprintf(stderr,"Unable to open file %s for writing\n",file);
              return EXIT_FAILURE;
         }
         for (int imom=0; imom<Nmom_eq; imom++) {
                for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust_real[iclust]=Zqr_eq[iclust][imass][imom][iOa2];
                fprintf(fout,"%8.10f %8.10f %8.10f  \n",SinP2_eq[imom],Zqr_eq[Nclust][imass][imom][iOa2],error_jack_real(aux_clust_real,Nclust));
          }
      fclose(fout);
     }



    free(Zq); 
	return 0;

}

